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Abstract 

A method to obtain (approximate) analytical expressions for the ra- 
dial distribution functions in a multicomponent mixture of additive hard 
spheres that was recently introduced is used to obtain the direct correla- 
tion functions and bridge functions in these systems. This method, which 
yields results practically equivalent to the Generalized Mean Spherical 
Approximation and includes thermodynamic consistency, is an alterna- 
tive to the usual integral equation approaches and requires as input only 
the contact values of the radial distribution functions and the isothermal 
compressibility. Calculations of the bridge functions for a binary mix- 
ture using the Boublfk-Mansoori-Carnahan-Starling-Leland equation of 
state are compared to parallel results obtained from the solution of the 
Percus-Yevick equation. We find that the conjecture recently proposed 
by Guzman and del Rio (1998, Molec. Phys., 95, 645) stating that the 
zeros of the bridge functions occur approximately at the same value of 
the shifted distance for all pairs of interactions is at odds with our results. 
Moreover, in the case of disparate sizes, even the Percus-Yevick bridge 
functions do not have this property. It is also found that the bridge func- 
tions are not necessarily non-positive. 
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1 Introduction 



Integral equation theories for the description of thermodynamic and structural 
properties of liquids usually lead to qualitatively satisfactory results. Nonethe- 
less, in general they involve hard numerical labor as well as criteria to formulate 
the closure relations that are not clearcut. This is true even for the simplest and 
most studied systems, namely the pure hard-sphere fluid and hard-sphere fluid 
mixtures. Therefore, it is not surprising that many attempts at providing gen- 
eral features, symmetries, approximations or parametrizations of the so-called 
bridge functions have been reported in the literature ||, [| U [5[ ^, fj], ||, ||, 
|Io| . These bridge functions enter in the closure relations and are defined 
as the sum of elementary diagrams (whose precise computation is a formidable 
and rather difficult task); they account for some molecular spatial correlations 
of higher order than pair correlations. Of course the availability of the ana- 
lytical results provided by the Percus-Yevick (PY) theory in the case of 
the pure hard-sphere fluid |l3| and hard-sphere fluid mixtures [|l4| allows one 
to determine explicitly the bridge functions in this instance, but they inherit 
the (theoretically) unpleasant lack of thermodynamic consistency as well as the 
limited density range of applicability involved in the PY approximation. 

In a related context, it is worth pointing out that in the pioneering work 
of Rosenfeld and Ashcroft |l5| it was found that an important class of pair 
potentials shared the property that their corresponding bridge functions were 
remarkably similar to each other and to the hard-sphere bridge function. This 
observation led to the common form of the reference-hypernetted chain theory 
[ fl6| , considered by many to be perhaps the most accurate theory for the struc- 
tural properties of fluids, in which the bridge functions of the system of interest 
are equated to those of a hard-sphere system. Thus, the search for accurate and 
relatively simple approximations for the bridge functions of a pure hard-sphere 
fluid and hard-sphere fluid mixtures has been pursued in the last few years. No- 
table among the results of this pursuit are the empirical parametrization due to 
Malijevsky and Labi'k (ML) for the hard-sphere fluid and its recent extension to 
binary hard-sphere mixtures [jl7| . These involve a careful and thorough analysis 
of a large set of computer simulation data. Some apparent regularities of the 
bridge functions in the case of binary mixtures (present in the PY results and 
in recent simulation data of Malijevsky et al. [Q) has been recently suggested 
by Guzman and del Rio ifjjfl . Were this regularities to hold in general, they 
would allow one to simplify the ML parametrization and serve as a starting 
point to consider hard-sphere mixtures with three or more components. Due 
to the scarcity of simulation data for mixtures, the suggestion remains a mere 
conjecture up to now. 

Notwithstanding the merits of all these theoretical and semiempirical efforts, 
it is clear that, especially in the case of mixtures, the scarcity of results and the 
relatively slow progress reflect the amount and difficulty of the numerical work 
that has to be done to get them. Therefore, one may reasonably wonder whether 
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an alternative theoretical approach, that at least avoided the inherent difficulty 
of solving nonlinear integral equations, would also provide information on the 
bridge functions of the pure hard-sphere fluid and hard-sphere mixtures. It is 
the major aim of this paper to provide an affirmative answer to the foregoing 
question. 

In previous work |2^, we have introduced a method to analytically derive 
(approximate) expressions for the radial distribution functions and structure 
factors of fluids and fluid mixtures. This method rests on a completely differ- 
ent philosophy than the one involved in integral equation theories and thus is 
totally void of the difficulty associated with providing any particular closure. 
In the case of pure hard-sphere fluids and hard-sphere mixtures and in the 
lowest order of approximation, it yields the well known PY results. Further- 
more and by construction, our expressions in the next order of approximation, 
which yields results practically equivalent to those of the Generalized Mean 
Spherical Approximation (GMSA) but is much simpler to implement, also 
embody thermodynamic consistency. As shown below, by using such radial dis- 
tribution functions in connection with the Ornstein-Zernike (OZ) equation it is 
rather straightforward to derive the direct correlation functions and, in turn, 
the bridge functions of the system, this latter for distances greater than the 
contact distance. 

The organization of the paper is the following. In Section ^| we outline 
the main ideas of our method to obtain the radial distribution functions of an 
iV-component mixture of additive hard-spheres. For this mixture the Boublfk- 
Mansoori-Carnahan-Starling-Leland (BMCSL) |23| equation of state and the 
Grundke-Henderson-Lee-Levesque (GHLL) |25[ contact values of the radial 
distribution functions, which yield the BMCSL equation of state, are assumed 
to hold. If N is set equal to one then the pure hard-sphere fluid case readily 
follows. Expressions for the direct correlation functions are given there. Section 
[| provides an analysis of the bridge functions as well as a comparison with 
previous work. We close the paper in Section |] with further discussion and 
some concluding remarks. 

2 The radial distribution functions, the direct 
correlation functions and the bridge functions 
of a hard-sphere mixture 

In this Section we outline the method to obtain (approximate) analytical expres- 
sions for the radial distribution functions (r) of a multicomponent hard-sphere 
mixture. It consists of an extension to mixtures of the method previously ap- 
plied to one-component systems of hard spheres, sticky hard spheres, and square 
wells [^0| . For details the reader may refer to Ref. |2lj . 

An iV-component mixture made of pi hard spheres (of diameter a{) per 



3 



volume unit may be characterized by 2N — 1 parameters (for instance, the 
N — 1 molar fractions Xi, the N — 1 size ratios Oijo\ and the packing fraction 
J2i Pi a f) an d involves N(N + l)/2 radial distribution functions gij(r). 
As happens in the PY and GMSA theories, it is convenient to work in the 
Laplace space and define 

/>oo 

Gijis) = / dre- sr r 9ll (r). (1) 
Jo 

There are two basic requirements that Gij (s) must satisfy. First, since gij (r) = 
for r < <Jij, with ov, = (cr^ + aj) /2, and gij{<J^j) = finite, this implies that (i) 
lim^oo s e SCTiJ 'dj -(s) = finite. Second, the isothermal compressibility kt = 
finite, so that (ii) lim s ^o[Gij(s) — s~ 2 ] = finite. The approximation we will use 
consists of assuming the following functional form: 

p — S<Tij 

Gij{8) = -2^^ Llk{s)[{1 + as)l ~ A ^]~V W 

k 

where I is the N x N unit matrix, 

M S )=^+LgWLgV (3) 

and 

2 

Aii(s) = Pi J2 Ms°iK +1 L%- n \ (4) 

n=0 

with 

<p n {x) = x-^+V 

Condition (i) is verified by construction. On the other hand, condition (ii) 
yields two linear sets of N 2 equations each, whose solution is straightforward, 
namely 

Lf) = A + \'<Tj + 2\'a - A Pk<TkL$ , (6) 

k 

L\V = + —OiOj + (A + \'ui)a - -a, ^ Pk^k^f] , (7) 

k 

where A = 2tt/(1 - rf) and A' = (A/2) 2 £ fe Pk°l- 

(2) 

The parameters L\j and a (which play a role similar to that of the param- 
eters Kij and z in the GMSA) are arbitrary, so that conditions (i) and (ii) are 

(2) 

satisfied regardless of their choice. In particular, if one chooses L\- = a = 0, 
our approximation coincides with the PY solution. If, on the other hand, we fix 
given values for gij{crfj), we get the relationship 

4 2) = 2naa ij9ij (a±). (8) 



E 

,m=0 



(5) 
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Thus, only a remains to be determined. Finally, if we fix kt, we obtain a closed 
equation for a of degree 2N . It is worth pointing out that in the particular case 
of a pure hard-sphere fluid (N = 1) one gets a quadratic algebraic equation for 
a, while for a binary mixture (JV = 2) the explicit knowledge of Gy(s) only 
requires to solve a quartic equation, which also has an analytical solution. A 
natural choice to close the scheme, which we will of course consider in this pa- 
per, is to take the GHLL values ||| of g i} ;(oy), as well as the corresponding 
BMCSL [^3) isothermal compressibility kj. But other possibilities are available 
and one of them will also be addressed later on. Once Gy (s) has been deter- 
mined, inverse Laplace transformation directly yields rgy(r), while the Fourier 
transforms hij (q) of the total correlation functions hij (r) readily follow from 
the relation 



hij(q) 



dr exp(iq ■ r)hij(r) 
2 _ Gij (s) - Gjjj-s) 



(9) 



s—zq 



where 1 is the imaginary unit. In Fourier space and introducing the quantities 
Hij( a ) — ^fWPjhij{q) an d Cij(q) — Jpipj Cij (?) the OZ equation reads 



C(g) = H( ? ).[l + H(flf)]- 1 , 



(10) 



so that after replacement of H (q) and subsequent inverse Fourier transformation 
it is straightforward to get Cij(r). The result gives Cy (r) for r > cry as the 
superposition of N Yukawas (see the Appendix): 



N R (i) 



C ij( r ) = eXP [ _2 ^( r ~ a 
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where q = ±izi with £ = 1, . . . , N are the zeros of det I + H (q) 
tudes K\;> are obtained by applying the residue theorem as 



(11) 

and the ampli- 



K m = iH- zl<Jii lim ZiM{q _ tz ^ 

•> Z7T q—>zze 



(12) 



Finally, we note that the bridge functions Bij(r) for r > Oij are linked to 
gij(r) and cy(r) through 



B i:j (r) = ]ng i:j (r) - g^r) + ClJ {r) + 1. 



(13) 



Equations ( |ll| ) and (|1^), after replacement of the results for gij(r), will be 
used below to investigate some properties of the bridge functions Bij (r) in pure 
hard-sphere fluids and hard-sphere mixtures. 
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3 Comparison with other results 



We begin with the pure hard-sphere fluid, that is, we now consider the case 
when N = 1. For this system a variety of closures to the OZ equation are 
available @, | g § g, §, for instance 

5 PY (r)=ln[l+ 7 (r)]- 7 (r), (14) 

B HNC (r)=0, (15) 

B VM (r) = [7(r)]2 (16) 

{> 2{l + a l7 (r)}' Ubj 

J B MS (r) = [l + 2 7 (r)] 1/2 - 7 (r)-l, (17) 

B RY (r) = In (l + ^P{[l-cxp(- a2 r)] 7 (r)}-l | _ 
[ 1 - cxp(-a 2 r) J 

B BPGG (r) = ^ + a37 ( r) ]Va3 _ 7(r) _ x ^ (19) 

where 7 (r) = #(r)-c(r)-l and the labels HNC, VM, MS, RY and BPGG denote 
the hypernetted-chain Verlet modified ||, Martynov-Sarkisov Rogers- 
Young Q and Ballone-Pastore-Galli-Gazzillo [|| closures, respectively, and the 
adjustable parameters ai(i — 1, 2, 3) have been estimated to take the values a 2 = 

0. 160 Q, ai = 4/5 and 03 = 15/8 0. It is at this point interesting to analyze to 
what extent the bridge functions calculated upon substitution in such closures 
of our expressions for the radial distribution function and direct correlation 
function differ or are compatible with the actual bridge function computed using 
equation (p"3]). To this end, we introduce the quantity AB*(r) = B*(r) — B(r), 
where the asterisk refers to a given label. Notice that in particular AB HNC (r) = 
—B(r) and AB PY (r) = ln[l — c(r)/g(r)\. In figure [l] we display the behavior of 
AB*(r) as a function of the shifted distance r — a for a packing fraction 77 = 0.49, 

1. e. close to the freezing transition. Clearly the main discrepancies between both 
types of calculation show up near the contact point, but one could reasonably 
argue that the RY, the MS and the VM closures are rather compatible with the 
result of the present approach in that region. 

Figure || displays the results for the bridge function as obtained with the 
PY theory, the parametrization of Malijevsky and Labfk (ML) jlTj and our 
formulation, again for the packing fraction r\ = 0.49. As clearly seen in the 
figure, the discrepancy between the results of the parametrization and both ours 
and those of the PY theory is rather significant. Also, although not perceptible 
in the scale of the figure, we note that in our case the first maximum of the 
bridge function attains a positive value, whereas both the ML parametrization 
and the PY theory always lead to non-positive values for the bridge function. 
We will come back to this point later on. 

We now turn to binary mixtures. In figures [3] and results for the different 
bridge functions are shown for two cases. In figure H, which corresponds to 
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the case examined by Guzman and del Rio |L9| , it is an equimolar mixture with 
77 = 0.35 and a%ja\ = 0.6; the second mixture is defined by 77 = 0.49, x\ = 1/16, 
X2 = 15/16 and o^/ci = 0.3. We have also included in these figures the results 
obtained with the PY approximation. The differences between the results of 
both approaches are particularly important in the first maximum of the bridge 
functions and once more in our case we get some intervals of positive values. 
These figures also illustrate the fact that the phase-shift symmetry that was 
recently conjectured by Guzman and del Rio Jl9| to hold on the basis of the 
behavior observed in figure [|, is not even present in the PY theory, particularly 
at short distances, and this becomes more evident as the disparity in size ratio 
is increased, as shown in figure [| 

As another illustration, in figures ^| and ^| we display results for a ternary 
mixture where 77 = 0.49, x\ — £2 = 1/102, 23 = 100/102, a%la\ = 0.3 and 
03/01 =0.1. Apart from exhibiting a more complicated structure than in the 
case of binary mixtures — notice, for example, the existence of a negative first 
maximum for B\\ and Byi — , again an important difference between our results 
and those of the PY theory is that we may get positive values for the bridge 
functions in some regions. 

Thus far we have only considered the GHLL prescription for the contact 
values giji&ij) and the isothermal compressibility kt derived from the BMCSL 
equation of state as the input in our method. One may reasonably wonder 
whether the use of different values for <7y(cr^) and/or kt would also yield sim- 
ilar results. In order to assess the importance of other reasonable choices, we 
have made calculations using the contact values gij (a^ ) obtained by extrapola- 
tion of simulation results |lS[ ] and the value of kt derived from an equation of 
state for mixtures (eCS) recently proposed by us ^6|. To illustrate the results 
one gets with these choices, in figure ^ we compare the various calculations of 
the function Bn(r) for the case considered earlier in figure^, namely the binary 
mixture defined by 7/ = 0.49, x% = 1/16, X2 — 15/16 and o^/o^ = 0.3. Except 
in the region up to the first maximum and near the second minimum, the curves 
obtained with either the BMCSL equation of state and with the second choice 
(eCS) using our procedure are practically indistinguishable. Nevertheless, al- 
though the first maximum with the eCS choice is still positive, its amplitude is 
much smaller than the one using the BMCSL equation of state, up to a point 
that the positive character can be hardly ascertained in the scale of the figure. 
It is also worth pointing out that the value of a is more sensitive to the choice of 
kt than to the values of ffij(c^). In fact, for this mixture one gets a = 0.0189 
using the BMCSL equation of state while a = 0.0118 using our proposal (eCS) 
for the equation of state. A smaller value of a means that cn(r) goes to zero 
more rapidly which in turn implies a much smaller (but still positive) value for 
the first maximum of Bn(r) in the latter case. 
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4 Discussion 

The points arising from the results of the previous sections deserve further elab- 
oration. To begin with, it is fair to say that our approach leads to a simple 
and clearcut procedure to determine both the bridge functions and the direct 
correlation functions in a multicomponent hard-sphere mixture, requiring as the 
only input the contact values of the radial distribution functions that specify 
the actual equation of state of the mixture. It should be pointed out that while 
the procedure is capable of yielding the values of the direct correlation functions 
for all distances, including those inside the hard cores, in the case of the bridge 
functions it is limited to distances greater than the contact distance. This is 
due to the fact that our method does not deal neither with closures nor with 
the cavity functions. Nevertheless, this restriction may be disposed of, at least 
for the case of the pure hard-sphere fluid, by considering approximate analytical 
forms of the cavity function that are available in the literature ||, . In this 

connection, we should mention that the form of the cavity function derived by 
Zhou and Stell || has been shown to be compatible with our g(r) in the sense 
that it yields the same values for both g{cr + ) and dg(r)/dr\ r=a+ |2^]. If the 
available simulation results for the radial distribution functions for mixtures are 
scarce, those for the bridge functions are to our knowledge nonexistent. In the 
absence of such data to compare with the results we have presented, it would 
be of course premature to reach definite conclusions. One could argue that the 
accuracy of the bridge functions might be estimated indirectly by comparing 
the radial distribution functions calculated using the OZ equation with a given 
closure and simulation results. In the present approach, however, this is unnec- 
essary since we have explicit (analytical) expressions for the radial distribution 
functions from the very beginning, and these have been already compared to 
simulation results both for the one-component and two-component hard-sphere 
systems in Refs. ^(J and pi]| . 

It is clear that a key difference between our results for the bridge functions 
and most of those previously reported, is the fact that in our case these functions 
may attain both positive and negative values. In connection with this issue 
one cannot overlook the fact that it has often been assumed that the bridge 
functions should be nonpositive. This is certainly the case in the PY theory and 
various parametrizations and approximations have included such an assumption. 
Nevertheless, as Rast et al. |ll[ have recently pointed out, there seems to be 
no rigorous reason or argument stating that it should be so. In fact, any theory 
that leads to a positive value of the direct correlation function c(r) at a distance 
where g{r) = 1 will produce a positive B(r) at that distance. For instance, 
taking the Monte Carlo data for c{r) obtained by Groot et al. [^9| and those of 
g(r) given by Barker and Henderson po| , one finds that B{r) ~ 0.2 for r ~ 1.85<r 
and -q = 0.445. In further support of the likely correctness of our results, one 
can invoke the fact that in our case thermodynamic consistency is an ingredient 
of the formulation while for instance the PY theory does not share this property. 
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For the pure hard-sphere fluid, the reasonable compatibility between our bridge 
function and the one computed using in particular the RY closure (which was 
originally proposed to achieve thermodynamic consistency) is very satisfactory 
in this respect. Further, the fact that our gij(r) are in better agreement with 
simulation results than those of the PY theory particularly in the region around 
the contact point also favors the present approach. 

Concerning the ML parametrization for mixtures and the apparent regularity 
of the shifted bridge functions that was conjectured to hold in general in Ref. 
[ fl9| , we can only add that unfortunately it does not do so. Indeed, it would have 
been rather remarkable that the relatively simple forms proposed for the bridge 
functions would have been able to capture the rich and varied behaviors that 
one would expect from the number of parameters involved in the description of 
mixtures. 

Finally, we want to point out that due to the similarity of the bridge functions 
corresponding to different potentials and those of hard-spheres — in fact the 
universality of the hard-sphere bridge functional has been recently shown to be 
very reliable JTo[ — these results are not only relevant for hard-sphere mixtures, 
but they may also prove useful in connection with the integral equation approach 
in liquid theory for mixtures with other interaction potentials. 

Fruitful discussions with F. del Rio, who suggested this problem to us, are 
gratefully acknowledged. One of us (M.L.H.) wants to thank the Universidad 
de Extremadura at Badajoz, where the draft of the paper was prepared, for its 
hospitality. He also acknowledges the partial support of DGAPA-UNAM under 
Project IN-117798. The work of two of us (S.B.Y. and A.S.) has been partially 
supported by DGES (Spain) through grant No. PB97-1501 and by the Junta 
de Extremadura-Fondo Social Europeo through grant No. IPR98C019. 



From the Fourier transform (q) one can get the direct correlation function in 
real space as 



i(r) 



1 



(2tt) 3 
1 



dq exp(-zq • r)cy(g) 
dqqe^Ctjiq). 



(20) 



It is now convenient to see q as a complex variable. Thus, if r > o~ij, it is 
possible to take a contour integration around the upper half plane in equation 
~20l) . According to Eq. (|l0|), the functions Cij{q) have a common set of poles, 



namely the zeros of D(q) = det I + H(g) 



A careful inspection of the results 



obtained from our method shows that the zeros of D(q) are the roots of a 



polynomial in q 2 of degree N. More specifically, the zeros of D{q) are q = ±izi, 
where the Zi {1 = 1,..., N) are positive real numbers. 
Application of the residue theorem yields 



N 

I 



c *j( r ) = n — y^^ e z ' r lim Cjj(?)(9-^)i r>aij, (21) 



=1 



which is equivalent to Eqs. ( p.l[) and (]12|). To be more explicit, let us rewrite 
Eq. © as 

C(g) = I - F- 1 ^), (22) 



where F(g) = I + H(q). Therefore, Eq. ( |i2| ) is equivalent to 

K? = ~^—= lim [F- 1 ^)] (g 2 + z, 2 ). (23) 
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Figure 1: The difference AB*(r) as a function of the shifted distance r — a for a 
simple hard-sphere fluid with rj = 0.49 and a = 1, according to various closures: 
HNC (dash-dotted line), PY (solid line), RY (dashed line), MS (dotted line), 
VM (solid line with triangles) and BPGG (solid line with circles). 
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Figure 2: Bridge function B(r) versus r — a for a simple hard-sphere fluid with 
r/ = 0.49 and a — 1. Solid line: present method; dashed line: PY results; dotted 
line: ML parametrization. 
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Figure 3: Bridge functions Bij (r) versus r—aij for an equimolar binary mixture 
of hard spheres with r/ = 0.35 and diameters er 2 = 0.6 and <7i = 1. Solid 
lines: present method; dashed lines: PY results. The curves for £?n(r) contain 
no symbols, those for i?i2(r) contain circles and the ones for B 2 2(r) contain 
triangles. 
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Figure 4: Bridge functions -By(r) versus r — Uij for a binary mixture of hard 
spheres with r] = 0.49, molar fraction x% = 1/16 and diameters 02 = 0.3 and 
o\ = 1 . The code for the different curves is as in figure 0. 
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Figure 5: Bridge functions Bij{r) versus r — for a ternary mixture of hard 
spheres with r\ — 0.49, molar fractions x\ = x 2 = 1/102 and diameters 03 = 0.1, 
(T2 = 0.3 and o\ = 1. Solid lines: present method; dashed lines: PY results. 
The curves for -Bn(r) contain no symbols, those for Bi 2 {r) contain circles and 
the ones for Bi 3 (r) contain triangles. 
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Figure 6: Bridge functions Bij{r) versus r — for a ternary mixture of hard 
spheres with r\ = 0.49, molar fractions x\ = x 2 = 1/102, and diameters 03 = 0.1, 
(72 = 0.3 and <ri = 1. Solid lines: present method; dashed lines: PY results. 
The curves for -622(7) contain no symbols, those for -623(f) contain circles and 
the ones for B 33 (r) contain triangles. 
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Figure 7: Bridge function Bn(r) versus r — o\ for a binary mixture of hard 
spheres with r\ = 0.49, molar fraction x\ = 1/16 and diameters 02 = 0.3 and 
(Tl = 1. Solid line: present method using the GHLL contact values of gij(cr^) 
and the compressibility kt derived from the BMCSL equation of state; dotted 
line: present method using the contact values of gij(ctj) obtained by extrapo- 
lation of the simulation data in Ref. Jlq| (namely, 511 = 10.23, gi2 — 4.69 and 
.922 = 3.57) and the compressibility kt derived from the equation of state eCS 
proposed in Ref. [^6); dashed line: PY results. 
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